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DIRECTION COSINE COMPUTATIONAL ERROR 


By John W. Jordan 
Electronics Research Center 


SUMMARY 

Strapdown inertial systems possess a potential advantage 
over more conventional inertial systems since they do not require 
a mechanical gimbal structure to maintain a stable coordinate 
reference. Instead, strapdown systems replace the gimbal struc- 
ture by a direction cosine matrix which is contained in and up- 
dated by the system computer. To assess the performance of a 
strapdown system or to determine the computer requirements , a 
detailed study of the computational error introduced by the com- 
puter is necessary. This report is concerned with the techniques 
which may be employed to estimate the computational error and its 
effects on system performance. Although focused upon a particu- 
lar problem of practical importance, many of the techniques have 
a wider applicability to aerospace software in general. The 
report is divided into the following sections: 

1. Introduction - The computer must solve the matrix differ- 
ential equation B = ftB where B is the computed direction 
cosine matrix and ft is a skew-symmetric angular rate matrix, 
the elements of which are determined by the system gyro- 
scopes. An error matrix (z) is defined which provides a link 
between the direction cosine error and system performance. 

The computer-generated error is divided into two categories: 

(a) Algorithm Error , which is due to the approximate nature 
of the numerical formulas used, and (b) Wordlength Error , 
which is due to a finite computer wordlength. 

2. Algorithm Error - A differential equation is derived for the 
propagation of the error (z) matrix. It applies to a general 
digital computational process with an arbitrary rate input. 
Since B should be an orthogonal matrix, it is possible to 
include a correction routine in the airborne computer program 
so that B will remain orthogonal despite computational error. 

A second differential equation is derived for the propagation 
of the z matrix when an orthogonality correction is used. 

Closed form analytical solutions are obtained for both 
error differential equations, provided that the angular rate 
matrix is self-commutative . 

A computable solution which allows the determination of 
the z matrix for any arbitrary angular rate input is developed. 
All solutions are verified by simulation. 



A compensation technique is devised (and verified by 
simulation) which greatly enhances the performance of low 
order algorithms. In many cases, this technique will reduce 
both the error and computer requirements . 

3 . Wordlength Error - Wordlength error is treated from the 
statistical viewpoint of Henrici. The theory is developed 
for linear matrix differential equations. The direction 
cosine problem is then used as a specific example. The re- 
sults of Henrici are extended to include different computer 
number systems (i.e., sign-magnitude or complement repre- 
sentations) and different organizations of the computer's 
arithmetic unit. The emphasis is upon a computable solution 
which may be applied to problems for which analytical solu- 
tions are either not possible or are not practical. The 
solution includes both fixed- and floating-point machines. 
The computable solution is verified by simulation. 

4. System Performance - The computable solutions provide the 
algorithm and wordlength error matrices for any arbitrary 
vehicle rate profile. It is possible to incorporate these 
solutions into a conventional error analysis program in 
order to determine overall system performance. The computa- 
tional error is handled in a manner completely analogous to 
the inertial instrument error and system performance is ex- 
pressed by identical criteria. 

A method (susceptible to automation) is presented for the 
"optimum" design of aerospace software. Although the approach 
is again general in nature, it is related to the problem of 
selecting the "best" direction cosine algorithm for a particular 
mission . 


I. INTRODUCTION 

The Apollo vehicle that lands on the moon will have as a 
backup system a " strapdown" inertial navigator. The term "strap- 
down" indicates that the accelerometers of the inertial system 
are rigidly connected ("strapped down") to the parent vehicle. 

The acceleration that they measure will be in body coordinates, 
i.e., dependent upon the instantaneous orientation of the vehicle, 
and not in the coordinate frame in which the navigation equations 
are solved. In order to remove the effect of vehicle rotations, 
a coordinate transformation matrix B(t) is calculated by the 
computer so that the measured acceleration vector in body coor- 
dinates a B may be resolved into the navigational frame as a^: 

a^(t) = B T (t)a B (t) (1) 
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The transposed matrix is used for convenience in the sequel. The 
matrix B(t) is 3x3 and contains nine elements called the direction 
cosines. Their time rate of change is given by 


B (t) = fi (t) B (t) 


(2) 


where £2(t) is a skew-symmetric matrix of body angular rates as 
measured by the system gyroscopes 


n (t) = 


o 

W y(t) 


w z (t) 

0 

-in (t) 
X 


-to ( t) 

y 

to ( t) 

X 

0 


(3) 


The computer solves the matrix differential equation (Eq. 

(2)) in real time. This computation must be done many times a 
second. To evaluate system performance, it is necessary to 
estimate the error introduced by the computational process. This 
report is not intended as a comparative analysis of the many 
procedures which have been proposed for solving the differential 
equations. Rather, it is an attempt to formulate a consistent 
analysis technique which may be applied to the different algorithms 
and missions. 

Error Criteria 

The direction cosines of a strapdown inertial system are 
used to solve the equation 

a N (t) = B T (t)a_ B (t) (4) 

where 

a B = the vector acceleration of the vehicle in body 
coordinates , 

a^ = the vector acceleration of the vehicle in naviga- 
tional coordinates, 

B = the direction cosine matrix which represents the 

transformation from body to navigational coordinates. 


At this point, the manner in which the computer solves the 
continuous Eq. (4) is not of concern. A digital computer will 
use a discrete formulation and this will introduce some error. 

To separate the effect of different error sources, it will be 
assumed that the computer is capable of solving Eq. (4) exactly, 
but that the direction cosines are not accurately known. In 
other words, the computer solves (exactly): 


» K - B T (t)a B (5) 

A 

where B(t) represents a "best estimate" of the true direction co- 
sine matrix B(t), and a N is the measured (not true) acceleration 
in the navigational frame. Defining 


E ( t ) = B(t) - B (t) , 


( 6 ) 


the acceleration error is given by 


«2„ = E a B 


(7) 


Since 


then 


-B B ^N ' 


6a_ T = E Ba._ 
— N — N 


Defining 


T 

Z = BE 


the two basic error equations are: 


e -b 

(8) 

T 

6<ir = Z a . 
— N — N 

(9) 


These equations provide two alternate descriptions for the effect 
of the direction cosine error. They are illustrated in Figure 1. 
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Figure 1.- Effect of direction cosine error 


Integration of the acceleration error described by either 
Eq. (8) or (9) for a specific mission profile will determine the 
errors in position and velocity which results from the direction 
cosine error. Although this is a fundamental technique for re- 
lating the direction cosine errors to system performance, it 
sometimes suffices to consider the cosine error matrices alone if 
one is not concerned with overall system analysis but simply the 
comparative performance of different cosine algorithms. 

As a matter of convenience, the Z matrix rather than the 
E matrix will be considered. Since any matrix is equal to the 
sum of a symmetric and a skew-symmetric matrix. 


Z 


z s + z ss 


( 10 ) 


with the symmetric part equal to 


- ,H z + zT ) 


(ii) 
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and the skew- symmetric part equal to 


Z 


SS 




(12) 


The Zgg matrix will contain three independent elements which will 
be denoted by 


V V Z Z • (13) 

The Zg matrix will contain six independent elements. The diagonal 
elements will be denoted by 


J XX' YY' ZZ 


(14) 


and the off-diagonal elements by 


Z 


XY' 


Z XZ' Z YZ * 


(15) 


Orthogonality Correction 

Since the direction cosine matrix B defines a transforma- 
tion of coordinates, it will be an orthogonal matrix. An ortho- 
gonal matrix has the property that the product of the matrix and 
its transpose is the unit matrix 


BB 


T 


I . 


(16) 


A 1 

The estimated cosine matrix B will not remain orthogonal because 
of computational error. Since Eq. (16) provides a simple check 
on the orthogonality of a matrix, it is possible to correct the 
non-orthogonality of B. The result will be a new estimate of the 
direction cosines B* which satisfies Eq. (16): 


Since an orthogonal matrix has the additional property of trans- 
forming a vector without changing its length, it would seem that 
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the new estimate would be superior to the old estimate for use in 
Eq. (5) since it will transform the acceleration vector without 
distorting its magnitude. 

The improvement (if any) in performance, because of the 
orthogonality correction, can be determined once the effect of 
the correction upon the error (z) matrix is established. 


The matrix B is contained in the computer as the best 
estimate of B. Since: 


B = B + E 
B = B (I + Z) . 


(16a) 


Let : 

T = ^B T B - i) ; (17) 

then 

T = Zg + J Z T Z . (18) 

From Eq. (17) it is obvious that the T matrix will be zero if B 
is orthogonal. If computational error renders the B matrix non- 
orthogonal, then the T matrix will be equal to the Zg matrix plus 
a second-order term. Error compensation may be accomplished by 
including a routine in the airborne computer which ensures that 
B remains orthogonal. Then the T matrix will be driven to zero: 


or 


1 T 

T = z s + I z 2 


= 0 


'sc 


1 T 

- i z T z 


(19) 


where the subscript SC denotes the value of the Zg matrix after 
compensation. Since Zg^ is a second-order term, 


Z 


SC 




( 20 ) 
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By substitution, it is easy to show that 

B* = B - BT (21) 


provides a new estimated cosine matrix which is orthogonal. 
Since the true cosine matrix is not known, approximate error 
compensation may be achieved by 


B* = B - BT (22) 


which differs from Eq. (21) by third-order terms. Non-orthogon- 
ality can be completely eliminated by an iterative application of 
Eq. (22) . 

Whether or not compensation is done depends upon the par- 
ticular mission requirements since a trade-off of software com- 
plexity versus accuracy is involved. 

However, it is clear that algorithms may be evaluated in 
either a compensated or uncompensated form. The skew-symmetric 
part of the Z matrix represents a more fundamental limitation 
since additional information is required to compensate for this 
error component . 

The use of an orthogonality correction will also affect the 
growth of the direction cosine error which arises during on-board 
solution of Eq. (2). This aspect of the problem has not been 
included in previous studies of the direction cosine computation 
and therefore is treated in some detail in the section on 
Algorithm Error. 


Computational Error 

The direction cosine matrix propagates as 


B (t) = ft (t) B (t) . 


Unfortunately, for any reasonable vehicle environment the elements 
of the ft (t) matrix will not be known a priori, but must be deter- 
mined in real time by monitoring the output of the system gyro- 
scopes. Therefore, an analytic solution is not possible in an 
operational environment and an estimated matrix B must be deter- 
mined by some computational technique. Current methods are based 
upon the use of standard numerical integration formulas to integrate 
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the matrix differential equations in real time. The computational 
error which arises during the integration may be divided into two 
categories : 

Algorithm error - which is due to the approximate nature 
of the numerical integration techniques used to integrate 
Eq. (2) . This error would exist even if the computer had 
infinite word length; 

Wordlength error - which is the additional error which 
arises from the finite computer word length. 

For many problems, these two error sources are essentially inde- 
pendent and may be evaluated separately. The total computational 
error will then be the sum of the two separate errors. Obviously, 
it would be a poor procedure to minimize one source of error 
without regard to the other. The first step in achieving a real- 
istic tradeoff between error sources, however, is an accurate 
assessment of the individual errors. 


II. ALGORITHM ERROR 

Algorithm error results from the approximate nature of the 
numerical integration techniques. The actual on-board integration 
is done in a discrete fashion. If the true direction cosine 
matrix propagates as 


b k+i - Vk <23) 

where is the state transition matrix between states K and K+l. 
Then the approximate cosine matrix propagates as 


B 


K+l 


$ B 
K K 


(24) 


Defining 


$ = $ - $ , 


subtraction gives 


E = 0 E + (ft R 
K+l K K V K K 


(25) 


9 


( 26 ) 


I 

for the stbp-by-step propagation of the error matrix. The Z 
matrix is i 


J K+1 


B K+1 E K+1 


K+l 


rn rn m m ^ m rn 

B $ $ B Z + B $ $ B Z + B $ * B 
K K K K K K K K K K K K K K 


While recognizing the fact that both B and $ are orthogonal and 
defining 


R* 


T T~ 

B $ $ B 
K K K K ' 


(27) 


then the equation for the Z matrix can be put in the form 


Z K+1 “ Z K = Az = r * z k + R * • (28) 

If h is the numerical integration step size (the time interval 
between steps) , then 

f - e( r * z k + R *) • (29 > 

The left-hand side now has the form of a difference quotient. 

The approximation 


Z 



(30) 


may be introduced. The term AZ may be regarded as the per-step 
error while 1/h is the number of integration steps per unit time, 
with the derivative being equal to the rate of growth of the Z 
matrix. Then 

Z = RZ + R (31) 

where 


R = R*/h = f c R* 


10 



and f c is the computational frequency in steps per second. Eq. 
(31) is the basic equation which describes the growth of the Z 
matrix because of the discretization error introduced by the 
digital computation process. The basic error equations are sum- 


marized in Table I. 

TABLE I 


PROPAGATION OF 
THE ERROR MATRIX 


b k+i $ k b k 


A A A 

b k+i = Vk 


A 

$ = $ - $ 


T T~ 

R = f B $ $B 
c 


Z = RZ + R 

Z (0) = 0 


Algorithm Error with an 
Orthogonality Correction 

The determination of the 
algorithm error is complicated by 
the use of an orthogonality cor- 
rection. The orthogonality of 
the estimated cosine matrix B can 
be tested at any time, and there- 
fore a correction for non-ortho- 
gonality can be made at any time. 

It is entirely feasible to update 
the estimated cosine matrix B 
for a number of cycles without 
using an orthogonality correction 
and to employ the correction 
routine only periodically to re- 
establish the orthogonality of 
the B matrix. To a certain extent, 
this procedure will correct the 
accumulation of past errors. 


The propagation of the Z matrix is given by 


Z = RZ + R . (31a) 


This equation may still be used when an orthogonality correction 
is included. Suppose that the correction is applied at the times 

t=nT n= 1,2,3,.... 


Then, at each one of these times the Z matrix must be re-initial- 
ized as 


Z*(nT) = Z ss (nT) 



(nT) Z (nT) 
2 



(32) 


The integration may then be continued to the next correction 
point . 
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For the greatest accuracy, the orthogonality correction can 
be applied everytime the estimated direction cosine matrix B is 
updated. This mode of operation, along with the completely un- 
compensated mode, provides upper and lower limits for the periodic 
use of an orthogonality correction. Therefore, it is desirable 
to determine the propagation of the Z matrix when an orthogonality 
correction is applied at every step. 

After the direction cosine matrix is updated, a correction 
matrix may be determined as shown previously 


T = 


b k+i b k+i 1 


or 


T = 


. (ilil) 


(i + 


,t~t; 

3 <J) 4 ^ 

K K K K 


K 


)- 


(33) 


Hereafter the subscript K will be suppressed. 

(i + z T )(i + 2r s + b t 5 t 5b)(i + z) - I 
T ~ ' 2 ~ 


(34) 


where 



(35) 


is the symmetric part of the R matrix. Neglecting the second- 
order term in Eq. (34) yields: 

(l + Z T )(l + 2R s )(l + z) - I 
T = 2 • (36) 


Interchanging the order of the factors (a second-order effect) 
gives 

[. I + 2Z C + Z T zUl + 2R ) - I 

T = V — — . (37) 

2 
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T 

One effect of the T matrix is to replace Zg by -Z Z/2 as noted 
previously. These two values may be considerably different if 
an orthogonality correction has not been applied recently. How- 
ever, if a correction is being used everytime the direction 
cosine matrix is updated, these matrices will be approximately 
equal in magnitude and 


2Z 2 + Z T Z = 0 (38) 

and 

T = R g . (39) 

i 

The propagation of the error matrix with a continuous orthogo- 
nality correction is therefore 


Z C = 


R ss z c + 


ss 


(40) 


where R^ s is the skew-symmetric component of the R matrix. This 
result is summarized in Table II. 


Closed-Form Solutions 


TABLE II 


Although the actual vehicle 

PROPAGATION OF THE ERROR MATRIX environment is too complex to 

WITH A CONTINUOUS allow a closed-form solution for 

ORTHOGONALITY CORRECTION the direction cosine matrix, it 

is a standard error analysis 
technique to assume an admittedly 
simplified model to make the 
analysis more tractable mathe- 
matically. Therefore, it will 
be assumed that the angular rate 
matrix n has a simple form which 
allows a closed-form solution 
for the direction cosine matrix 
(B) to be found. Appendix A shows that it is always possible to 
take the initial condition of the B matrix as the unit matrix 
without loss of generality. This initial condition will be 
implicit in the analysis that follows. The reader is reminded 
that the significance of the E and Z matrices is due to Eqs . (8) 

and (9) -- that is, they provide a link between the cosine errors 
and overall system error. For the time being however, we are 
only concerned with determining the functional form of the E and 
Z matrices. 
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There is one class of closed-form solutions which is easily 
found. It is well known that the system of matrix equations 


X = A ( t ) X 


X(0) = I 


has the unique solution 


s: 


X (t) = exp I A(X) dX , 


provided that 


A(t 1 )A(t 2 ) = A(t 2 )A(t 1 ) 


(41) 


(42) 


(43) 


for all tj and t 2 . Reference 1 proves that this equality is 
satisfied if, and only if, the matrix A(t) can be represented by 


A (t) = Z i K i f i (t) (44) 

where the fj_ are scalar functions and the K-^ are mutually com- 
mutative constant matrices. For skew-symmetric matrices this 
condition is satisfied when the constant matrices are proportional. 
Therefore, the angular rate matrix has the form 


fi(t) = Kf(t) (45) 

where 

f(t) = Zf i (t) (46) 


The statement is self-commutative" will be used to indicate 
that Eq. (43) is satisfied for all t and that the angular rate 
matrix has the form of Eq. (45). Then, the equation 

B = QB B (0 ) = I (47) 
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has the solution 


Next, if 


J“ 


B = exp K | f (X ) dX 
'0 


( 48 ) 


and if 


/‘ 


F (t) = | f(X) dX 

: 0 


B = e 


KF 


( 49 ) 


K = 


0 

-k 


, k 

L y 


o 

-k 


-k 


x 


and 

oj n — k + k + k~ 

0 x y 2 

2 

so that ojq is the sum of the squares of the elements in the 
constant matrix K, then 


K 3 = 


K 4 = 


K 


2 v 2 

-ojq K 


( 50 ) 


Using these identities all higher powers of the K matrix may be 
reduced to either the first or second power. When the exponential 
function is expanded and like powers of the K matrix are collected, 
the solution for the B matrix assumes the form 
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B(t) = I + 


sin w^Fft) 


K + 


(51) 


U3 


0 


1 - cos a)QF(t) 2 
2 K 


This is a closed-form solution for the cosine matrix under the 
restriction of proportional angular rates about each coordinate 
axis. To calculate the algorithm error it is necessary to find 
the appropriate value of the R matrix in Table I. The true state 
transition matrix for one computational step from t-h to t is 
given by 


and if 


$ (t) = exp (KF(t) - KF(t-h)} , 


H ( t) = F (t) - F (t-h) , 


then 


sin u.H 1 - cos w n H „ 

<S(t) = I + — K + * — K 


OJ, 


Cjj, 


(52) 


(53) 


The R matrix becomes 

R = f c «£i K (54) 

because of the commutivity of the matrices. The state transition 
matrix $ is given by Eq. (53) , while the approximate matrix $ 
depends on the numerical integration algorithm which is being 
evaluated . 

Since all powers of the K matrix can be reduced to the 
first or second power, it follows that the general form of the 
R matrix is 


R 


r x K + r 2 K 


2 


(55) 


where r^ and r 2 are scalar time functions. 

2 

Obviously, the K and K matrices commute, and using the 
theorem quoted from reference 1, the corresponding solution for 
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the Z matrix is 


S i K SoK 2 

Z = e e -I 


(56) 


where 



dX 


dX 


(57) 


(58) 


By expanding Eq. (56) and collecting terms, the general form of 
the Z matrix 


Z = Z ss + Z s (59) 

can be written as 

Z = z ss K + z g K 2 (60) 


where z gs and z g are the scalar functions 


2 

-a) 0 S 2 

e sxn Oq s^ 


J SS 


0) 


0 


Z S = 


n -a) 0 S 2 

1 - £ COS COq S^ 


a ), 


(61) 


The error relationships are summarized in Table III for the 
special case of proportional input rates. For this case the nine 
separate elements of the Z matrix need not be calculated, since 
the two scalar parameters z g and z completely describe the 
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TABLE III 


PROPORTIONAL RATES 



the error. As an example the solution for 


f2(t) = K^a^ + a 2 t + a^ sin w s tj 

is easily found since all the elements of the Q matrix are pro- 
portional. When this restriction is met, an analytical solution 
for the Z matrix requires only the integration of the R matrix. 

When a continuous orthogonality correction is used, the 
symmetric part of the R matrix is compensated for. The solution 
for the Z matrix may be modified by setting r 2 (and therefore 8) 
equal to zero. The results are summarized in Table IV. Note 
that the compensated solution differs from the uncompensated 
solution by the absence of the exponential term. Whether or not 
this is of practical importance depends upon the magnitude of the 
exponential term. 


Numerical Integration Formulas 

A large number of numerical integration formulas have been 
proposed for the integration of direction cosines. Since this 
report is concerned with analysis techniques rather than with a 
comparative analysis of the different formulas, one or two of 
the currently used algorithms will be selected for use in examples. 


At this point, it will be assumed that the inertial system 
is instrumented with rate-integrating gyroscopes . These gyroscopes 
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TABLE IV 


PROPORTIONAL RATES 
CONTINUOUS ORTHOGONALITY CORRECTION 



have an output which is porportional to the integral of the input 
rate or 


0 (t) 



fl(X) dA 


(62) 


where 0 (t) is a skew symmetric matrix of gyro outputs. 
When 


fl(t) = Kf (t) , (63) 

then 

0 (t) = KH (t) (64) 

and 

$ = e 0 . (65) 


Approximate state transition matrices may be generated by ex- 
panding the exponential function in a Taylor series and retaining 
only the higher order terms 

e 0 = I + 9 + 0 2 /2 + 9 3 /6 + ... (66) 
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First-Order Taylor Series: 


$ = I + 6 (67) 


Second-Order Taylor Series : 


8 = i + e + e 2 /2 


( 68 ) 


Note that these approximations do not assume that the angular 
rate matrix A is constant during an integration step. The basic 
assumption is that A is self-commutative. There is, of course, 
the additional judgment that the number of terms retained is 
sufficiently accurate for the application at hand. 

Narrow-band Solution 

Although it is possible to determine an exact form of the 
R matrix for many representative rate matrices, it is much simpler 
to use an approximate solution which is adequate for almost any 
practical situation. For this reason, exact solutions will be 
discussed in Appendix B, while the remainder of this section will 
consider an approximate solution. 

Given the state transition matrix 


sin a> n H 1 - cos oo n H _ 

0 = 1 + — K + = — K 


03 


(69) 


0 


co 


0 


where 


H = F ( t ) - F(t-h) , 


a narrow-band solution may be defined for the situation where 


a) Q H < 0.5 . 

With this restriction 

sin WqH = WqH - 


(70) 
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COS 0) Q H = 1 


+ 


2„2 

W 0 H 


4 4 
W 0 H 
24 


When the equality of Eq. (70) is satisfied, the magnitude of the 
first neglected term is only one tenth that of the fourth-order 
term which is retained. The term "narrow-band" is applied since 
the functions sin oigH and cos cogH also arise in the study of 
narrow-band FM signals . Then 


$ = I + 



(ji 


2 IT 
0 6 



oi. 


H_ 

24, 


K 


(71) 


It is now a simple matter to determine on approximate R matrix 
for the Taylor Series algorithms. 

First-Order Taylor Series 

S=I+0=I+ HK (72) 

R = f 

c 

R = -(*c“o h3/6 ) K - ( £ c h2/2 ) k2 

Second-Order Taylor Series 

2 2 

$ = I + 8 + = I + HK + K (73) 

R = ( f c oigH 3 / 6) K - (f c ai^H 2 /s)K 2 . 


From the mean value theorem of differential calculus: 


H ( t) = F(t) - F(t-h) = hf(t-r)h) 


where 


0 < n ^ 1 • 
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Therefore : 


H (t) = hf (t - rjh) . 


The mean value of f (t) may be approximated by f(t) itself since 
the small shift in the time axis of the error solution is not 
significant. The final form of the R matrices for the Taylor 
Series algorithm are listed in Table V. 


TABLE V 

TAYLOR SERIES NARROWBAND SOLUTION 


1st Order 

(-oo 2 /3f 2 )f 3 (t)K - (l/2f c )f 2 (t)K 2 

2nd Order 

(co 2 /6f 2 )f 3 (t)K - (co 2 /8f 3 )f 4 (t)K 2 


Table V then provides a value of the R matrix for any arbitrary 
rate input, provided that the rate matrix ( t) is self-commutative, 
and provided that the narrow-band restriction 


u) Q hf (t) < 0.5 


is satisfied. 

Note that the narrow-band criteria are not dependent upon 
the frequency content of the rate matrix, but rather upon its 
amplitude. For example, the rate inputs 


co = co = co = — cos t (74) 

X y Z ^ 


CO = CO = CO = — cos lot 

x y Z /3 


both have the same amplitude with 


co 


0 


10 , 
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although their frequency content differs by a factor of 10. If 
the computational rate is 20 times a second or greater 


oj q H <0.5 


and both inputs may be categorized as narrow-band. For the 
20 sample/sec situation, the inaccuracy in the error calculation 
would be approximately 10 percent. However, for any application 
with such high amplitude input rates the sampling frequency would, 
in all probability, exceed 20 samples/sec. This would reduce the 
inaccuracy in the error calculation. The narrow-band criteria 
will adequately represent almost any practical situation. 

Example I — Constant Rate 

Let 


Then 


f(t) = 1 
SHt) = K . 

H (t) = h . 


From Table V, the R matrices are 
First-Order : 


R “ -(“0 /3f o) K - ( 1/2f c) K2 (75) 

Second-Order : 

R = ( W 2 /6£ 2 )k - (^/8f 3 )K 2 . (76) 

In both cases the R matrix is a constant and the integration 
called for in Table V is easily performed. The results appear as 
Table VI . 



TABLE VI 


CONSTANT RATE NARROW-BAND SOLUTION 
TAYLOR SERIES ALGORITHMS 



Example II — Sinusoidal Rate 
Let 


f (t) = cos to t . 

s 


Then 


2 /hco 

H ( t ) = — sin ( — =— | cos 
w \ z 

s \ 


[-.(■ - 5)1 


( 77 ) 


which for small h is approximately 


H (t) = h cos w s ^ / 


( 78 ) 


so that 


^ = w s 2 


which may be approximated by 


H (t) = hcos a) t . 

s 


The R matrices are then 
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First-Order 


R = -(w 2 /3f 2 ) cos 3 w s t K - (l/2f c ) cos 2 to g t K 2 ( 

\ 

Second-Order 

R = (wg/6f 2 ) cos 3 <u t K - (tD 2 /8f 3 ) cos 4 w g t K 2 ( 


These expressions for the R matric'es must now be integrated in 
accordance with Table III to obtain the error terms z gs and z g 
The results are given in Table VII. 


TABLE VII 

SINUSOIDAL RATE NARROWBARD SOLUTION 
TAYLOR SERIES ALGORITHMS 


1st Order 


fa/3f 2 W * ln ^ 



6 

2nd Order 

a 




For both the sinusoidal and constant rate cases the total error 
is 


8 sin a ,, , 1 - e p cos a T .2 


E 








When a continuous orthogonality correction is used, the error 
becomes 


Z 

c 


sin a 


K + 


1 - cos a 



K 


( 82 ) 


Simulation 

A solution for the Z matrix can often be found by simula- 
tion. The algorithm under study is programmed on a digital 
computer and used to integrate a specified 0 (t) for which an 
exact cosine matrix is known. The exact solution is then sub- 
tracted from the simulation result to determine the error matrix 
(see Figure 2) . In this report the term "simulation" will always 
be used in this context. Its purpose will be to provide a check 
on the analytical results . 



Figure 2.- Effect of direction cosine error 

Simulations were conducted for both sinusoidal and constant 
rate inputs. The results are shown in Figures 3 through 6. 

Parameter values for the simu- 
TABLE VIII lations are given in Table VIII. 

The error was calculated for 40 
SIMULATION PARAMETER VALUES or 400 seconds. The figures 

compare the simulation results 
to a plot of the analytical 
solution obtained using the 
narrow-band approximation. 

These cases are repeated with a 
continuous orthogonality correc- 
tion in Figures 7 through 10. 


03 = 

X 

0.3 

f = 32 

03 = 

y 

0.3 

c 

to = 0.25 

03 = 

Z 

0.3 

s 
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SIMULATION RESULT ANALYTICAL SOLUTION 


Figure 3.- Constant rates, first-order Taylor series 
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10 100 200 300 400 


10 IOO 200 300 400 


SIMULATION RESULT ANALYTICAL SOLUTION 

Figure 7.- Constant rates, first-order Taylor series/ 
continuous orthogonality correction 






SIMULATION RESULT 


ANALYTICAL SOLUTION 


Figure 8.- Constant rates, second-order Taylor .series , 
continuous orthogonality correction 
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The difference between the analytical and simulation results 
of Figure 9 is of interest. In deriving the differential equation 
(Table II) for the growth of the Z matrix with an orthogonal cor- 
rection, it was assumed that 

2Z + Z T Z = 0 ; (83) 

s 

that is, with the correction continuously applied the B matrix 
will remain orthogonal. The compensated Z matrix is then 


Z 

c 


z K 
ssc 


z K 
sc 


(84) 


for which Eq. (83) is satisfied. However, this approximation 
will be in error by the non-orthogonality error introduced since 
the last orthogonality correction. More accurately 


Z 

c 


z K 
ssc 


+ (z 


sc 


z )K" 

p 


(85) 


where Zp is the error introduced since the last correction. 

When the orthogonality correction is applied Zp will be compensa- 
ted. However, it will generate a second-order residual because 
of the Z T Z term. This residual error results from the incomplete 
compensation of the symmetric component of the Z matrix by the 
orthogonality correction. When an orthogonality correction is 
applied at each integration interval, Zp will be equal to the 
symmetric error introduced in a single step. From Table V it is 
apparent that while this error may be appreciable for a first- 
order algorithm, it will be much smaller for the second-order 
algorithm. The residual error is completely insignificant in the 
simulation results of Figure 10 for the second order algorithm. 

The residual error is not serious because it is not propa- 
gated. Each time the orthogonality correction is applied it 
compensates for the old error and introduces a new residual 
error. The new error will depend on the amplitude of the angular 
rates over the integration interval. The residual error can be 
reduced for the first-order algorithm by decreasing the integra- 
tion interval, and thereby reducing the per step error. It can 
be eliminated by applying the orthogonality correction twice per 
integration step. The second application corrects for the 
residual error of the first. Figure 9 contains a plot of the 
symmetric error term when an orthogonality correction is used 
twice each integration step to correct for the residual error. 
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z ssc 



z ssc 




SIMULATION RESULT 


SIMULATION RESULT 


Figure 11.- Constant rates, first-order Taylor series, 
periodic orthogonality correction 



For the angular rates, integration interval and simulation 
times are chosen, and the parameters a and 3 remain small. The 
solutions of Tables II and III may be approximated by 


z 


ss 


— (1 
“0 


+ 3) 


( 86 ) 


z 


ssc 


a 



z 

s 


2 


a 


2<al 


■d + 3) 


z 


sc 



The difference between the compensated and uncompensated 
solutions is due to the parameter 3 which, in turn, is propor- 
tional to the symmetric component of the R matrix. The first- 
order algorithm is characterized by a large symmetric error 
component. It is therefore useful for short-term missions which 
involve low angular rates. The effect of 3 on the accuracy of a 
first-order algorithm is evident when a comparison is made between 
Figures 3 and 7 . 

A more interesting comparison might be made between a com- 
pensated first-order algorithm and an uncompensated second-order 
algorithm. At the same computational frequency both should re- 
quire roughly equal computer time and memory. A sinusoidal rate 
about one or more axes can be used to approximate vehicle limit 
cycling or angular vibrations. As illustrated by Figures 6 and 9, 
the first-order algorithm would have a skew-symmetric error twice 
that of the second-order algorithm. This error is oscillatory. 

The uncompensated second-order algorithm, on the other hand, has 
a secular symmetric error term, while the symmetric error in the 
first-order algorithm is bounded. The choice of a "best" algorithm 
is not always simple since it involves a trade-off between the 
allowable error, the expected environment, and the computational 
resources available. A subsequent section will consider this 
trade-off in more detail. 

In Figures 7 through 10, the orthogonality correction was 
applied at each integration step or 32 times a second. It is 
possible to apply the correction much less frequently. In 
Figure 11 an orthogonality correction is periodically employed. 


37 


A first-order algorithm with a constant input rate was used. 

The correction was applied every 6 4 integration steps (once every 
2 seconds) . Since the residual error of a first order algorithm 
would be quite large after 64 steps, the correction was applied 
twice at the 2-second intervals. By comparing Figures 3 and 11 
it is evident that even the infrequent use of the correction pre- 
vented the exponential increase in the skew-symmetric error term. 

Figure 11 actually provides a plot of the symmetric error 
immediately after a correction has been applied. (The error is 
plotted at 10-second intervals over a total time of 400 seconds.) 
Additional plots in Figure 11 show the error for the first 10 
seconds at intervals of 0.25 second. Note that the error grows 
in an uncompensated manner between the 2-second compensation 
points. Therefore: 

1. Figure 3 represents the uncompensated error, 

2. Figure 7 represents the continuously compensated error, 

3. The 400-second plot of Figure 11 is somewhat misleading since 
it represents the symmetric error immediately after an ortho- 
gonality correction. This discrepancy may be resolved by 
also plotting the error immediately before the orthogonality 
correction. The effective error may then be considered to be 
the average of the "before" and "after" values. The effective 
error is then directly dependent upon the frequency at which 
the correction is applied. This frequency is another trade- 
off parameter between computer loading and error level. 

General Solution 

When the angular rate matrix is self-commutative, i.e.: 


ft(t-^)ft(t2) = (t2) (t-^) 


(87) 


for all t and R matrix is immediately known (Table V) . In addi- 
tion, if ft (t) can be expressed as the sum of simple analytical 
functions, it is usually possible to integrate the R matrix and 
find the total analytical solution for the algorithm error 
(Table III). In general, however, ft(t) will have the form: 


ft(t) = 2K i f i 


( 88 ) 


where the Kj matrices are not mutually commutative. In this case, 
a complete analysis is much more difficult, since no analytical 
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solution for the direction cosines is known. The well known 
Neumann series does provide a symbolic solution, however: 


/' 


B (t) =l+| fl(X) dX + / I dxdX + 


"0 


V 


( 89 ) 


Therefore : 


//. 


$(t) = 1+1 Si(X) dX + I I n(X)fi(T) drdX + ... (90) 

t-h *'t-h*'t-h 


from Eq. (62) 


/: 


0 (t) = | fl(X) dX . 

-h 


Now define: 


r 


0 n ( t ) = | fi(X)0 n _ 1 (X) dX e 0 (t) = I ; (91) 

t-h 


then 


$ = 


I + 0 + 



dX 


This series may be written as 


0 


i + e + 


if 

n=2 *t-h V. 


de 


n 


dt 


+ 0 


n-1 


1 

n! 



(92) 
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or 


$ 


oo 



n=2 


( 93 ) 


where 


C n 



nl 


d8 n \ 
dt J 


dA 


(94) 


Of course, Eq. (93) is no simpler than the original Eq. (89). 
Actually, it is more complex. However, Eq. (93) provides a link 
with previous error estimates. If the angular rate matrix is 
self-commutative, the right-hand side of Eq. (93) reduces to the 
exponential function alone. The series in C* may thus be inter- 
preted as the additional error resulting from a non-commutative 
rate matrix. The exponential function may be written as 


= I + 


sin 4) (t) 
<j> (t) 


0(t) + 


1 - COS (j) (t) 

<J > 2 (t) 


e 2 (t) 


(95) 


2 

where <f> (t) is a scalar function equal to the sum of the squares 
of the three independent elements of the matrix 0 (t) . Eq. (95) 
is a generalized form of Eq. (51) which represented the exponen- 
tial function for the case of self-commutative rate matrix. The 
error due to the approximation of the exponential term may be 
determined by the same procedure which was applied to the earlier 
situation. Since 


let 


rn rn 

R = f B r$B , 
c 


w = f c $ T $ = f (£ - ei) 


(96) 
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Once again 


0 (t) = hft (t - nh) 
0 < n < l 


by the mean value theorem. To simplify the determination of the 
narrow-band approximation may be used. The result is presented 
in Table IX: 


TABLE IX 

TAYLOR SERIES NARROWBAND SOLUTION 
GENERAL RATE MATRIX 


Order 

W 

NB 

First 

-w 2 (t)ft(t) , 7 

2 2 F~ (t) " C 2 (t) _ C 3 (t) 

3 f c 

c 

Second 

00? (t) ft (t) oo 2 ft 2 ( t) 

— 5 ^ 0 C,(t) - C-(t) 

6f 8f 

c c 


Table IX is a generalization of Table V for the case where ft (t) 
is not necessarily self-commutative . The two matrices C2 and C3 
represent the additional error which results from the non- 
commutativity of the angular rate matrix. This error is commonly 
called the "commutatively" error. From Eq. ( 96 ) : 


= f C* 
r c 2 


C 3 “ f o (C 3 - SC 2 > 
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From Eq. (94) : 


C* = 
z 


if 

J t- 

r 

J ±.1 


fi(A)0 (A) - 0 (A)J1(A) dA 


(97) 


C* = I n(A)0(A) - [e 2 (A)fi(A) + 0(A)fi(A)6(A) 
t-h 

+ (A) 0 ( A ) /6 J dA 


(98) 


Since 0 and are skew-symmetric, the principal commutativity 
term will be also. All three matrices have three independent 
elements. Eq. (97) may also be written in terms of a vector 
cross-product 


* 



1 

2 



0x(jO dA 


(99) 


where each vector contains the three independent elements of the 
corresponding matrix. This equation allows another interpretation 
of the major commutativity error term; it results when the vectors 
0 and a) are not co-linear during the integration interval. This 
implies that the direction of the angular rate vector is changing 
with time. 

Although a non-commutative rate matrix considerably compli- 
cates an error analysis, it is, unfortunately, the rule rather 
than the exception. The solution for a self-commutative rate 
matrix is very useful for a comparative analysis of different 
direction cosine algorithms; however, when a more realistic 
environment is assumed to determine system performance, non- 
commutativity errors inevitably arise. In these circumstances, 
it is necessary to evaluate their contribution to the total error. 

Example I 


The trajectory of a launch vehicle may often be described 
by a series of constant rate turns. It follows that an analytical 
solution is known for the direction cosine matrix. Since the 
turning rate is generally small, almost any numerical integration 
algorithm will demonstrate very small errors. However, a more 
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realistic system performance analysis might include the addition 
of sinusoidal rates about one, or more, axis to represent angular 
vibrations or vehicle limit cycles. For simplicity, consider 
only one sinusoidal term: 


ft(t) = Kq + cos to g t (100) 

Most likely, the matrices Kq and will not commute and a com- 
mutativity error will result. Let: 


L = KjKq - (101) 

The commutativity term of most interest is C 2 since it is of 
highest order and being skew-symmetric will not be reduced by an 
orthogonality correction. A direct evaluation of Eq. (97) gives: 


C 2 (t> 


c 

= ~ L(P 3 cos 0 J s t + P 2 sin to s t) 


( 102 ) 


where 


2 sin to h h(l + cos to h) 

P = S s 

2 2 


(103) 


to 


to. 


P 3 = 


h cos to s h 2(1 - cos w g h) 


to. 


to 


If: 


to h « 1 
s 


the C 2 matrix is approximately 


to L 

C 9 (t) = =- sin to t . (104) 

2r S 

c 
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E, 


In this case, the commutativity error term is sinusoidal as 
might be expected from the interaction of a constant rate and a 
sinusoidal rate. 

Example II 


The so-called "generalized coning motion" is defined by an 
angular rate matrix, the elements of which are 


to (t) 
X 

= to „ 

xO 


toy(t) 

= to . COS to t 

yO s 

(105) 

to z (t) 

= to A sin to t. 
zO s 



The principal commutativity term is most easily found from 
Eq. (99). If the error vector is denoted by 


then 


f 2 = 


(106) 


(i) „(j) „ 

C = _£0_JE0 f p 
x 2 c 1 


(107) 


tO „to ,, 

xO zO c (r , . , . . . 

c = f (P„ cos to t + P- sin to t) 

y o c 2 s 3 s 


w „to_ 


= — — f (-P„ sin to t + P_ cos to t) 
2 c 2 s 3 s 


where 


, sin to h 

P _ J2 s_ 

1 to 2 

s to 


( 108 ) 


and P 2 and P^ are defined in Example I 
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Again, when 03 s h « 1 


c = 


03 ..oi _03 

yO zO s 
2 

12f 

c 


03 -0) -03 

xO zO s . 

c — COS 03 t 

y i2f 2 s 


03 -.03 -03 

c = — Y — — — s i n (jj t 

z 12f 2 3 

c 


( 109 ) 


The x axis commutativity term is a constant. A constant term in 
the R matrix will cause a secular term in the Z matrix (see 
Table I) . This secular term may greatly degrade overall system 
performance. The commutativity terms for the y and z axes, being 
periodic, are much less important. 

Example III 


To determine a general form for the commutativity error, 
note that, by the mean value theorem: 


where 


fl (A) = (t Q ) + fi (t 1 ) A 


t Q = t-h < t ± < t , 


0 < A < h 


( 110 ) 


although the theorem does not provide the value of t^. 
Let 


= fi(t Q ) (111) 

= Q, (t^) . 
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A straightforward evaluation of Eq. (97) gives 


C 2 “ Q °^ l) 


(112) 


and 


c* = — K (ft,^ - n^n,) 

24f 


C 2 f c C 2 12f 2 ^ 1 S 0 fi Q^l ) 


(113) 


C 3 24f 3 ( ^i Q 0 2 ^o^l fi O + fl 0^1* 


It follows that both C 2 and C 3 are skew-symmetric. Since 
C 3 will usually be much smaller than C 3 , it may be neglected in 
comparison to C 2 . For both first- and second-order Taylor series 
algorithms, the commutativity error may be considered to be skew- 
symmetric and inversely proportional to f 2 . 

Once the commutativity error matrix is found, the W matrix 
is easily obtained from Table IX. Unfortunately, even if the W 
matrix is easily integrable, a closed-form solution for the Z 
matrix is not available as it was in the case of a self-commutative 
rate matrix. However, a solution may always be computed by numeri- 
cal integration techniques . A procedure for determining the Z 
matrix once the W matrix is known is outlined in the section on 
Performance Analysis. It will not be pursued further because a 
subsequent section will demonstrate that a computable solution for 
the Z matrix which does not require the determination of the W 
matrix is possible. 


Compensation Technique 
Given the basic error equation [Eq. (31)]: 


Z = RZ + R 
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the R matrix can be expressed as the sum of a symmetric and a 
skew-symmetric matrix: 


R = R + R . (114) 

ss s 

For a short time, or while Z is small, the initial error growth 
can be approximated by 


Z = R 


(115) 


if an orthogonality correction is not used, and by 


Z 

ssc 


R 


ss 


(116) 


Z = R Z 
sc ss ssc 


when the correction is used. Continuing with the case of no 
orthogonality correction, let 


S 



= S 


ss 


+ S . 


s 


(117) 


If B a is calculated using a second-order algorithm, the 
initial error growth will be 


B = B (I + S 
A ss 


+ V 


(118) 


Sg S is inversely proportional to f£ and S s 


where 

O O O 

proportional to 
algorithm but with 
calculate B A . 


is inversely 
Assume that Bg is calculated with the same 
computational rate one half that used to 


Then 


B d = B ( I + 4S + 8S ) 
B ss s 


(119) 
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These two relationships can now be used to eliminate either S ss 
or S s . 

For example : 


B 

c 



( 120 ) 


The skew- symmetric term is compensated at the cost of an increase 
in the symmetric term. This term, however, is often much smaller 
than the skew-symmetric term. A corresponding correction for a 
first-order algorithm is 


B 

c 



/I + 2S 


B l- 


( 121 ) 


which provides a reduction in the symmetric term as well. 

Of course, it is possible to use an orthogonality correction 
in the calculation of B A and Bg. However, the orthogonality cor- 
rection may be more economically applied to B c instead of to B A 
and Bg individually. Since this forces 


S = S 2 /2 , (122) 

s ss / ' 


it follows that S s will be compensated for, as well. If this 
procedure is periodically repeated, both the skew-symmetric and 
symmetric error will be compensated. The long-term error propa- 
gation will still be given by 


Z 

c 


R Z + R 
SS c ss 


(123) 


only R ss will now be an uncompensated higher order term. For a 
second-order algorithm, this term will be inversely proportional 
to f J and the algorithm will have a fourth-order error character- 
istic. A first-order algorithm will have a higher order symmetric 
term inversely proportional to f^. Since this term can be elimi- 
nated by an orthonormality correction, a fourth-order error 
characteristic can also be achieved. The commutativity error, 
being inversely proportional to f2, will also be compensated for. 
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Hopefully, the computer memory requirements are not con- 
siderably increased if the same routine is used to calculate B* 
and Bb* It is not necessary to store Bq since it is an updated 
B^. In addition, the time required to compute Bb will be only 
half that required to compute B^ because of the lower repetition 
rate. Since the error level is now that of a fourth-order 
algorithm, it is possible to reduce the basic computational 
frequency. For many applications the result will be an increase 
in accuracy, a decrease in the overall computational time, and a 
minor increase in memory requirements. 

The compensation technique is illustrated in Figure 12. 

Note that the correction frequencies may be submultiples of the 
computation frequency. Therefore, there are many variations 
possible, with each yielding a different accuracy versus computer 
loading trade-off. 

Figure 13 is the result of a simulation using the compensa- 
tion technique described in this section. The simulation used a 
second-order Taylor series. One run was made without an ortho- 
gonality correction and one was made with an orthogonality cor- 
rection applied to B c at each step. The matrix Bb was also set 
equal to B c after each correction. Figure 13 may be compared to 
Figures 6 and 10 since they represent the same angular rate matrix 
and computational frequency. 

A constant rate case is shown in Figure 14. This run may 
be compared to Figure 3 and 7 since they all represent a first- 
order algorithm with the same input and computational frequency. 

The first-order compensated solution of Figure 14 may also 
be compared with the second-order results of Figure 8. The 
second-order solution (using a continuous orthogonality correction) 
had a skew symmetric error of 


z - 0.2x10 
ssc 


-1 


The completely compensated first order method of Figure 14 had a 
corresponding error of 


z =0.9x10 
ssc 


Table X contains a rough estimate of the relative computer 
time requires for both algorithms. The total time is approximately 
the same for both. The fully compensated algorithm, however, has 
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NO ORTHOGONALITY CORRECTION 


Figure 13.- Sinusoidal rates, 
compensated solut: 
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NO ORTHOGONALITY CORRECTION CONTINUOUS ORTHOGONALITY 

CORRECTION 



CORRECTED RESIDUAL ERROR 


Figure 14.- Constant rates, first-order Taylor series, 
compensated solution 
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TABLE X 


■!' 

1 


"V-i. 


COMPUTATION TIMES 


Calculation 

2nd Order 

1st 

(compensated) 

b a 

1 

0.5 

b b 


0.25 

B c 


0.25 

Orthogonality 

Correction 

1 

1 

Total 

2 

2 


a much smaller error. If the computation frequency of the first- 
order algorithm were decreased by half, its error would be in- 
creased by a factor of 16 resulting in 


z = 0.14x10 * . 
ssc 

The compensated algorithm allows a 50 percent reduction in the 
computational load and still retains an order-of-magnitude 
accuracy improvement. 

Computable Solution 

The fact that B is an orthogonal matrix may be used to allow 
the separation of the error terms in and B 0 . 

Let 

T2 = B^B b - I . (124) 

For a second-order algorithm without an orthogonality correction: 
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B A = B(I + S ss 


(125) 


+ V 


b B = B(I + 4S ss + 8S s } 


T2 = 3S + 9S + (S e - S ) (8S o + 4S ) 
ss s s ss s ss 


T2 = 3S + 9S , 
ss s 


provided the second order terms are small 
Then 


or 


T2 - T2 


T 


ss 


„ _ T2 + T2 ' 

s " i 8 


and for a first-order algorithm 


T2 - T2 ' 


ss 


S = 
s 


T2 + T2 


T 


When an orthogonality correction is used, the higher 
terms cannot be neglected. The value of T2 becomes 


T2 = 3S + 17S 
ss s 


4S 


2 

ss 


(126) 

(127) 

(128) 

order 

(129) 


where the third term on the right-hand side is the most signifi- 
cant of the higher order terms. It is still true that 



ssc 


T2 


(130) 


S 


- T2 


T 


6 


however, the calculation of the symmetric error term is more 
complex. For a second-order algorithm 


and 


sc 


S 2 /2 
ss / 


sc 


T2 + T2 
18 


(131) 


(132) 


Equation (132) also applies for a first-order algorithm, provided 
the residual error is small. If the residual error is large, 
then 


and 


S > S /2 
sc ss 


sc 


T2 + T2 ' 
34 


(133) 


is a better estimate of the symmetric error. The factor of 18 is 
suitable for error analysis purposes since it will give a pes- 
simistic estimate (by a factor of 2) for a first-order algorithm 
when the residual error is large. It is also possible to use 


T2 - T2 ' 


ssc 


S 

sc 


T2 + T2 T , „2 / 8 \ 
32 ss\ 32/ 


(134) 


which is more complex, but has the advantage o,f being applicable 
to both first- and second-order algorithms, regardless of the 
residual error. 
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Since S ss and S s approximate Z ss and Z s , these equations 
may be used to find a computable solution for Z ss and Z s . The 
solution employs a FORTRAN program. The algorithm under study is 
a subroutine in this program. An arbitrary rate input is inte- 
grated to obtain Ba* The matrix Bg is computed at half the com- 
putational frequency of Ba/ using the same subroutine. From Ba 
and Bg the matrix T2 is found and then Z ss and Z s . 

Unlike the case of the compensated solution (discussed in 
the previous section), the Z matrix is not corrected. Therefore, 
the S and Z matrices will differ as time passes. The inaccuracy 
in the computable solution can be estimated from the Numman series 
solution for the Z matrix 


Z (t) 


S (t) 


+ 



R(X) S (X) 


+ . . . 


(135) 


The second term can be used to estimate the inaccuracy or to 
apply a correction. A simpler technique for checking the validity 
of the computable solution is to recompute the error with a 
different step size and check for a proportional increase in the 
error. 


A simulation program was written using the classical coning 
motion: 


y 

II 

o 

• 

3 




w 

s 

= 0. 
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w 

X 

= w 

s 

(1 - 

cos y) 
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= w 
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y 

cos 

w t 
s 

w 
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= w 
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sin 

Y 

sin 

w t 
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This is the only case (known to the author) where an analytic 
direction cosine matrix is known for a non-commutative rate matrix. 
A solution is contained in References 2 and 3. A summary of the 
work reported in reference 3 is contained in reference 4. The 
errors determined by simulation were compared to the errors deter- 
mined by a computable solution. The results are shown in Figures 
15 through 18. Since the rate matrix is not self-commutative, 
the error cannot be adequately described by the two scalar para- 
meters z ss and z s . The nine separate parameters defined in the 
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I 10 20 30 40 


COMPUTABLE SOLUTION 


-order Taylor series, continuous 


Error Criteria section must be used. For simplicity, only three 
of the nine are drawn in the figures. 


III. WORDLENGTH ERROR 
Introduction 

One of the most important factors influencing the design or 
evaluation of an airborne computer is an appropriate choice of 
wordlength. Although many factors must be considered in selecting 
the wordlength, the necessity of achieving adequate computational 
accuracy is often of decisive importance. Unfortunately, it is 
often quite difficult to determine the degradation of computational 
accuracy because of finite wordlength. 

Simulation is a useful tool for this purpose. A large, 
general-purpose computer can be programmed to simulate the air- 
borne computer and to execute mathematical calculations with a 
reduced wordlength. However, this approach is expensive, not 
only because of the computer running time, but more importantly, 
because of the programming effort involved. The basic limitation 
of simulation, however, is the difficulty in interpreting the 
simulation data without the benefit of an analytical theory. The 
wordlength error generated in a sequence of calculations is a 
function of the variables involved and will depend on the initial 
conditions. Changes in the initial conditions often generate 
very different error time histories. In these cases, it is a 
statistical description of the error which is of value. These 
statistics can be estimated by taking an ensemble average over 
many experiments, but, needless to say, the simulation costs may 
be prohibitively increased. 

The development of an analytical theory likewise presents 
many difficulties. The wordlength error depends upon: 

1. The order in which calculations are preformed, 

2. The operation of the computer's arithmetic unit, 

3. The fixed-point scaling used or the use of a floating 
point , 

4. The way negative numbers are represented in the 
computer , 

5. The initial conditions of the problem. 
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6. The way numbers are shortened to the required wordlength, 

7. The use of double precision or extended precision opera- 
tions . 

Two different approaches have been taken in generating an 
analytical theory. One approach is exemplified by Wilkenson's 
book (ref. 5) on determining maximum bounds for the error. An- 
other approach is used by Henrici (ref. 6). It deals with the 
statistics of the errors. In a recent paper Hull and Swenson 
(ref. 7) present the results of simulations designed to test the 
validity of the statistical model. They conclude that the models 
are, in general, very good. Discrepancies are both rare and mild. 
The discrepancies were attributed to the capability of the com- 
puter to do some operations exactly, whereas the statistical model 
assigns an error to all operations. Therefore, the results of the 
statistical model tend to be conservative. Others (refs. 2, 6, 

8, 9, 10) have also reported good correlation between the theory 
and simulation results. In the sequel, the statistical theory of 
Henrici (ref. 6) will be emphasized. To simplify the discussion, 
the following sections will consider the propagation of wordlength 
error for a set of first-order differential equations. The results 
will then be applied to the direction cosine problem. 

Propagation of Wordlength Error 

Many navigational and guidance functions can be described 
by a set of first order differential equations 


x = f (x, t) 


(137) 


The differential equations will be assumed to be linear and 
of the form 


x (t) = F (t) x (t) 


(138) 


where F(t) is, in general, a time-varying matrix. For convenience 
in digital computation, the differential equations can be con- 
verted into linear difference equations of the form 




= *k2k-i 


( 139 ) 


x^. is the n dimensional state vector, and 
is the nxn state transition matrix. 
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The digital computation will not, however, generate the 
true state vector x k because of the error introduced by the finite 
computer wordlength. Instead, a computed state vector x* will 
result 


x* = x* , + e. x* = 0 (140) 

— k k— k-1 — k — o 

where 

x£ is the computed state 

£ k is the n dimensional wordlength error introduced at 
the kth stage 

Now, let 



be the state error vector. Subtracting Eq. (139) from (140) 
gives : 


— k = »kSk-l + £-k (142) 

for the propagation of the state error. The statistical theory 
of wordlength error assumes that is a random vector, the value 
of which at any time is unknown. The mean value and variance of 
the error are assumed to be known. The mean value of the state 
error is 


e, = $. e, t + e, 
— k k— k-1 — k 


(143 


Let : 

m k = -e k (143) 

a, = 

— k — k 
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Then 


E k - 'A-1 + 2 k 


rn = 0 
o 


(145) 


Eq. (145) is the basic equation for the propagation of the mean 
value of the state error. To determine variations about the 
mean value, let 


— k ~ -k -k 


(146) 


2k " ^k " 2 k 


It follows from Eq. (143) that 


2k = $ k2k-l + 2 k - 


(147) 


and the covariance matrix of the variation about the mean value 
is 


p k * HjA 


(148) 


Writing 


P k ■ ®k P k-l*k + V^k + Vk-l^k + <Vk®k 


Q “ ‘Sjc'lk 


and making the assumption that 


Vk-A + %Vi'k = 0 ' 


(149) 


the covariance matrix becomes 


P. = 4> P. , Qp + Q. ; P =0 
k k k-1 k k o 


(150) 
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Equations (145) and (150) are two key relations since they allow 
a recursive calculation of the mean and covariance of the state 
error resulting from the wordlength error. They involve only the 
statistics of the error and do not require simulations to deter- 
mine the actual errors. 


Continuous Model 

Although the computer calculates in a discrete fashion, it 
is also possible to propagate the wordlength error by a continuous 
model (ref. 6). In this case the basic difference equation 


= \£k-i 


(151) 


is replaced by the differential equation 


x(t) = F ( t ) x ( t ) 


(152) 


where F(t) is in general a time-varying matrix. Since the 
relationship between the discrete and continuous model is well 
known (refs. 10,11,12) the key equations are listed in Table XI 
without derivation. 


TABLE XI 
KEY EQUATIONS 


Discrete Model : 

— k " \Ek-l + 2k m o “ 

p k - Vk-l'I + Q k P o 

0 

= 0 


Continuous Model: 



m(t) = F(t)m(t) + f c a(t) 

m (0 ) = 0 


P(t) = F (t) P ( t) + P(t)F T 

(t) + f^Q(t) 

P(0) = 0 
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In the continuous model of the mean value equation, the 
mean value of the wordlength error is given by f c a(t). The 
parameter f c is the computational frequency with which the dis- 
crete equations are being solved by the computer, i.e., the 
number of wordlength errors generated per unit time. 

Both the continuous and discrete models provide a computable 
solution for the error statistics rather than a closed form solu- 
tion. For some applications, the continuous model can be used to 
find closed-form solutions. Table XII contains the symbolic 
solution for the continuous model. 


TABLE XII 

SOLUTION FOR THE CONTINUOUS MODEL 



The results contained in Tables XI and XII are essentially 
those given by Henrici in reference 5. To use either the con- 
tinuous or discrete model, it is necessary to obtain expressions 
for the mean value and covariance of the word length error vector. 
In the next section the theory contained in reference 6 will be 
extended to include different computer number systems (i.e., 
complement representations as well as sign-magnitude) and different 
organizations of the computer's arithmetic unit. The extension 
will first be applied to fixed-point computers. A separate section 
will then extend the results to floating-point machines. 
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Fixed Point Computer 



Origin of wordlength error .- Wordlength error may result 
whenever the number of bits representing a number is reduced. 

Each computer operation must be examined individually to deter- 
mine whether or not it involves a reduction in significant bits. 
First, however, the process of shortening a single word will be 
analyzed . 

Positive numbers .- Consider a positive fixed point number x 
which is t bits in length. It is to be shortened to £ bits. The 
process of shortening may be done in two different ways, each of 
which generates errors in different ranges. 

Chopping . - The t bit word is simply chopped off to £ bits. 


4 t » 


f/AM. 91 

•* £ * 

m d ► 


In the chopping operation, d bits are lost which have some 
positive value r. If x lies in the range 0 < x < 1, then the 
maximum value of the error introduced by chopping will be 


Max r = -2 


In general, x will be scaled at some value 2 V , so that the maxi- 
mum value of the error is actually 


_v-£ 

Max r = -2 


The minimum value of the error "zero" occurs when the lost bits 
have the value "0". The actual value of the error lies some- 
where between these two limiting values. The statistical theory 
of wordlength error assumes that the error is uniformly distribu- 
ted between the extremes. Letting k = 2 V “ ^ , the probability 
distribution of the error is shown in Figure 19. For positive x 
the wordlength error e is simply the negative of the value of 
the lost bits or as shown in Figure 19. 
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Figure 19.- Distribution 
of wordlength error when 
a chopping operation is 
performed on a positive 
number . 


The error distribution shown in Figure 19 has an average value 
of 



This characteristic of the error caused by chopping is often 
highly undesirable. The second method of shortening wordlength 
does not have this drawback. 

Symmetrical rounding .- Before the t bit word is shortened 
to SL bits , the value k/2 is added. The wordlength is then short 
ened to SL bits. 

The rounding operation shifts the distribution of the word 
length error to the right so that it has zero mean value. 


* 


p(€) 





Figure 20.- Distribution of wordlength error when 
a symmetrical rounding operation is used. 
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Negative numbers .- For the case where x is a negative num- 
ber the results are slightly more complex. The sign of the word- 
length error will depend on the method of representing negative 
numbers in the computer, i.e., sign and magnitude, one's comple- 
ment or two's complement representations. (Positive numbers have 
identical representations for all numbers systems.) 

Sign and magnitude .- In this system, a negative number is 
represented by its absolute value plus a sign bit. A chopping 
operation will always reduce the absolute value of the number. 
Therefore, if the original number x is negative, the chopping 
operation will result in a less negative number or a positive 
error. Thus, the error for both positive and negative x can be 
written as 


e = -r sgn(x) 


(153) 


where 


sgn(x) =0 x = 0 

sgn(x) = -1 x < 0 


sgn(x) = 1 


x > 0 . 


For negative values of x the probability distribution of the 
error shown in Figure 19 is reflected about the vertical axis 
since the error will always have a sign opposite to the sign of 
x. A one's complement machine will have the same error character- 
istics as a sign and magnitude machine. 

Two's complement.- According to reference 13, the signed 
two's complement representation of a negative binary number will 
give the true value of the binary number if the value of the sign 
bit is regarded as negative. For example, the number 1,1001 
has the value 


-10000 + 1001 = -0111 


A wordlength error will result in a more negative result 
giving a negative error. Therefore, when negative numbers are 
represented in two's complement form, the wordlength error caused 
by chopping will always appear as in Figure 19 and will be 
independent of the sign of x. 


69 



Symmetrical rounding .- Finally, since symmetrical rounding 
produces an error distribution which is symmetrical about zero, 
the sign of x is of no consequence irrespective of the number 
system employed. 

Statistics of the wordlength error .- The variance of a 
variable uniformally distributed between 0 and k is given by 



The mean value and variance of the wordlength error are summari- 
zed in Table XIII. Use is made of the definitions 

a - -e and q = e 2 


TABLE XIII 
ERROR STATISTICS 


Number 

System 

Chopping 

Symmetric 

Rounding 


a 

q 

a 

q 

S/M 

j sgn(x) 

k 2 

0 

k 2 

2C 

k 

2 

12 

u 


Computer mechanization .- The errors introduced by the 
fundamental arithmetic operations may be analyzed in terms of 
the shortening of wordlength discussed in the previous sections. 

Although the system is described by the difference equation 


*k = $ k*k-l ' (154) 

this is not the usual form of implementation on a digital com- 
puter. A more common procedure involves the numerical integration 
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of the basic differential equations. This process may be written 
as 


2k - 2 k _i + *2 k (155) 

where the increment Ax is determined by the numerical integra- 
tion formula. To preserve accuracy, the step size is chosen small 
enough so that Ax is a relatively small quantity. Five different 
methods of performing the addition are of interest: 

SPSA - Single Precision Short Accumulator 
SPLA - Single Precision Long Accumulator 
DPSA - Double Precision Short Accumulator 
DP LA - Double Precision Long Accumulator 
EP - Extended Precision 


Single precision short accumulator. - Since Ax is small 
relative to x, assume that it is calculated with a smaller scale 
factor: 

2 V scale factor of x 
2^ scale factor of Ax 
s = v - y . 


Then Ax must be shifted s bits to 


to x. 


k-1 


the right before it is added 


x 


k-1 



x 


k 



The number Axjc is shortened to 5,-s bits. For sign and 
magnitude representation, the maximum error would be 

max e = -2^ ^ s sgn (Ax) = -2 V ^ sgn (Ax) . 
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The maximum value of the error depends not on the scaling of Ax, 
but rather on the scaling of x. The sign is opposite that of Ax. 

It has been assumed that the computation of Ax was carried 
out with an optimum scale factor. Wordlength error associated 
with the computation of Ax would then affect only the least sig- 
nificant bits of Ax and would be "shifted off" when the addition 
to x is performed. It is possible to program so that an explicit 
shifting of Ax is not required. In this case, the shift operation 
is performed by a multiplication during the calculation of Ax. 

Single precision long accumulator .- Assume that all vari- 
ables are of single precision with a length of A bits. However, 
the accumulator into which x^-^ is loaded is extended with 
additional bits which are equal to zero 


A k-1 

Ax 


The final result x^ is shortened by chopping or symmetrical 
rounding. For sign and magnitude arithmetic the maximum error 
would be 

max e = -2 sgn(x) . 

The magnitude of the error is the same as for single pre- 
cision addition, but the sign is now determined by x rather than 
Ax. For a two's complement number system, or if symmetrical 
rounding is used, long accumulator addition does not differ from 
short accumulator addition. 

Extended precision addition .- For extended precision addi- 
tion, the variable x is stored in the computer by two A bit words. 
Even if only the most significant A bits are used in other calcu- 
lations, the double length word is available when Ax is added to 
x k-l : 


A k-1 

Ax 
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The wordlength error results from shortening Ax and for a sign 
and magnitude number system is 


max e = -e^ x sgn(Ax) 


where e^ x is the wordlength error generated in calculating Ax. 

It must be estimated for each separate application. The magnitude 
of the error depends on the scaling of Ax as does its sign. 

It is not necessary to think of extended precision opera- 
tions as requiring a double length word to represent x. One may 
consider x as a single length word. Extended precision addition 
then requires a separate memory location to accumulate the lower 
bits of Ax which would be "shifted off" in doing a single-precision 
addition to x. It is then a matter of personal preference whether 
or not the memory location which accumulates the lower order bits 
of Ax are considered to be a double length extension of the vari- 
able x. 

Double precision addition .- If all the variables are repre- 
sented in double precision form (2£ bits) the word length error 
will be the same as for single precision addition with £ replaced 
by 2£ . 


Summary . - The error caused by the different forms of addi- 
tion are summarized in Table XIV. For a two's complement number 
system, the errors differ, at most, by a constant. In the table, 
it has been assumed that the error in calculating Ax is bounded 
by 


'Ax 


< 2 


w 


The value of w is highly program-dependent and must be estimated 
for each individual application. 


Floating Point Computers 

In a floating point computer, the scaling of the variables 
is itself a variable. The scaling of the variable x is given by 


2 


v— 1 


< X < 


2 


v 


where v is a computed quantity. The error statistics of Tables 
XIII and XIV are valid for floating-point as well as fixed-point 
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SPSA/DPSA 


Number System 
k 

a (chop) 

a (round) 

2 

q 


SM 

Two ' s 

Complement 

2 V “ £ /2 v ~ 2£ 

2 v-l /2 v-2Z 

k/2 sgn(Ax) 

k/2 

0 

0 

k 2 /12 

k 2 /12 


SPLA/DPLA 


Number System 

SM 

Two ' s 

Complement 

k 

2 V “V2 V ~ 2£ 

2 v-£^ 2 v-2 ^ 

a (chop) 

k/2 sgn(x) 

k/2 

a (round) 

0 

0 

2 

q 

k 2 /12 

k 2 /12 

EP 

Number System 

SM 

Two ' s 

Complement 

k 

2 w-£ 

2 w-£ 

a (chop) 

k/2 sgn ( x) 

k/2 

a (round) 

0 

0 

2 

q 

k 2 /12 

k 2 /12 











computations with the exception that v and w must be regarded as 
variables for the floating point calculations. If the equations 
of Table XI for the propagation of the wordlength error are solved 
on a computer, the fact that v is a variable represents only a 
slight complication. For example, on a digital computer, it is 
only necessary to include a subroutine which accepts a state 
variable as its input and determines 2 V for that variable. 


Closed-Form Solutions 

Covariance of the error .- Closed-form solutions can often 
be found by performing the integrations indicated in Table XII. 

The solution of the covariance equation is the same for all 
number systems. It is the complete solution if symmetric rounding 
is used since the mean value of the error will be zero. A 
simplification results if the initial reference time t 0 is equal 
to zero, since the first term in the formula for the covariance 
vanishes . 

Mean value of the error .- The mean value need only be 
computed when a chopping operation is used, since it is zero when 
a rounding operation is employed. Complications are introduced 
for a sign and magnitude number system since the sign of the 
error will be a variable. The integration must be divided into 
time periods which have the same error sign and a piecewise inte- 
gration carried out. 


Computable Solutions 

The value of closed-form solutions lies in the insight they 
provide into the propagation of the wordlength error. For all 
but the simplest systems, however, the closed-form solutions can 
involve tedious algebraic manipulations . These can be avoided if 
a computable rather than a closed-form solution is employed for 
analysis . The computable solution may be solved on a general- 
purpose digital computer or on an analog computer. 

Discrete computable solution .- A generalized computer pro- 
gram may be written which propagates the difference equations: 


x. = x, . 
— k k— k-1 


= *k^k-l + 2k 


m Q = 0 


P. = $. P. , + Q. 

k k k-1 k k 


P Q = 0 


(156) 
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The initial state Xq and the state transition matrix are 
inputs to the program. If $ varies with time, then a subroutine 
may be added to calculate the state transition matrix. The 
quantities and are calculated in accordance with Table XIV. 
Some of the quantities in Table XIII depend on the current state 
xj.. For example, if sign and magnitude numbers are assumed, the 
wordlength error will depend on sgn(Ax). In this case, the sub- 
routine calculating will use two values of the state vector 
to calculate: 


= " 


— k-1 


(157) 


The reader is cautioned that the function sgn(x) is defined with 


sgn(O) = 0 


The sign function found in most FORTRAN system is similar to sgn 
but it may have the value ±1 when its argument is zero. It will 
then be necessary to test for a zero argument before using the 
sign function. 

When a floating-point number system is being analyzed, the 
value of k in Table XIII will be a variable. This presents no 
problem for the computable solution. The current value of k may 
be found from the current state by a masking operation. A mask 
is applied to the variable x^ so that its exponent is left un- 
changed, while its mantissa becomes the number +1. The masking 
is done by a FORTRAN logical operation and no assembly language 
programming is necessary. The output of the program is a plot 
of X]^, m^, and P. The value of P may be interpreted as the 
variance of the error which would result if a rounding operation 
were used. If chopping is used, the error is given by the plot 
of the mean value of the state error (m) while P represents the 
variance of variations about the mean. 

Continuous computable solution .- A continuous solution may 
be found by solving the equations 


x(t) = F ( t) x (t) 
m ( t ) = F(t)m(t) + f c a(t) 


P(t) 


F (t) P ( t) + P ( t ) F 1 (t) 


+ f^Q(t) 


( 158 ) 
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on a digital or analog computer. This form has the advantage 
that the state transition matrix need not be known. It allows a 
solution for cases where no analytical state transition matrix 
can be found. In addition, the continuous model may be inte- 
grated on a digital computer with a variable step size algorithm. 

Summary . - Either the discrete or continuous computable 
solutions offer a very general and simple method for determining 
the propagation of the wordlength error. Table XII shows that 
the mean value and standard deviation are both linear functions 
of f c and 2 _ &. Therefore, the dependence on the parameters f c 
and % is known. In addition, it is usually possible to obtain 
some analytical results for at least a simplified model of the 
system of interest. These factors usually ensure that a high 
degree of confidence can be placed in the computable solution. 

Direction Cosine Wordlength Error 
The propagation of the direction cosine matrix is given by 


B = QB . 


The fact that a digital computer has a finite wordlength will 
result in the calculation of an inaccurate cosine matrix denoted 
by B. Although the error is random in nature, it is possible to 
calculate the mean value and variance of the error. 

Error propagation with symmetric rounding .- If the computer 
employs a symmetric rounding operation, the mean value of the 
state error vector will be zero and need not be calculated. The 
actual error will have a random nature which is characterized by 
the covariance matrix P(t) of the state error vector. Since the 
B matrix consists of direction cosines , the maximum value of any 
element is one. (Actually, the variables must be scaled slightly 
larger to allow for computational error; however, unity scaling 
will suffice for an error analysis.) 

Let the state variables be the first column of the B matrix 
and assume that a fixed-point computer is being used. If the 
wordlength is SL bits, then from Table XIV 


k = 2 


Since all state variables have the same scaling, the Q matrix 
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(Table XIII) is: 


Q = 



( 159 ) 


where I is the unit matrix. This value applies to single pre- 
cision and double precision (with Z replaced by 21 ) . With 
extended precision addition the magnitude of the error is pro- 
portional to the scaling of the state variable increment instead 
of the scaling of the state variable. The value of k may then 
differ by a power of 2. The general form of the solution will be 
the same. A closed-form solution for the covariance matrix can 
be found from Table XII . 


P(0) 
P (t) 

P (t) 


f 2 | B (x) Q (x) B T (x) dx 


0 
I 


f 2 k 2 

- ™ j B(x)B T (x) dx , 


and, since B is an orthogonal matrix: 


P(t) = 


2 2 
f k 
c 

12 


tl 


C160) 


(161) 


The covariance matrix of the error vector, because of finite 
wordlength, remains diagonal and grows linearly with time. The 
errors in the state variables remain uncorrelated. A similar 
analysis can be made for the two remaining columns of the B 
matrix. Because of the identical scaling, the three columns of 
the error matrix [the E matrix of Eg. (6)] will be identical. 
Since the error is a statistical quantity, it is necessary to 
generalize the definition of the error matrix. 
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Defining : 


U 



1 

1 

1 


1 

1 

1 


( 162 ) 


the E matrix may be defined as 


E = crU 


where 


a = -2- Vt 

yi2 


( 163 ) 


is the standard deviation (one sigma value) of the error. (Two 
or three sigma values can also be used.) Because the errors are 
uncorrelated, the Z matrix is also 


Z = crU . 


When symmetric rounding is used, the resultant errors are un- 
correlated with a standard deviation that grows with the square 
root of time . 

Mean value of the error .- If the computer wordlength is 
shortened by a chopping operation , the mean value of the error 
will, in general, be non-zero. The propagation of the estimated 
direction cosine matrix is given by 


B = QB + A 


( 164 ) 


where A is a matrix containing the average value of the wordlength 
errors introduced in one computational step. From Table XI, it 
follows that 


E = ftE + A 


( 165 ) 
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and 


• T 
Z = BA 


describe the propagation of the mean value of the direction 
cosine errors. Define the sgn function of a matrix as 


sgn ($) = [sgn (<j> ± j)] 


( 166 ) 


so that sgn ($) is a matrix derived from the matrix $ by replac- 
ing each element of $ by the values -1, 0, +1, depending on 
whether the element is less than, equal to, or greater than zero, 
respectively. Then Table XV contains the A matrix for several 
different number systems. The value of k is given for a fixed- 
point computer. When a floating point computer is used, k 
becomes a variable. 


TABLE XV 

WORDLENGTH ERROR MATRIX 


Number 

System 

S/M 

SPSA 

S/M 

SPLA 

2 1 s Comp 

U = 

‘l 1 1“| 
111 
_1 1 1J 

A 

^ sgn ($) 

^ sgn(A$) 


CM 

II 

44 


Closed-form solutions.- Closed-form or approximate solu- 
tions for the mean value of the wordlength error are possible. 

For example, the mean value of the two's complement wordlength 
error is 

Z T = U I B (x) dx . (167) 

*'o 
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For a constant rate about each axis : 



fl = K 


( 168 ) 


and 


B = 


sin w n t (1 - cos w n t) 0 

I + r JL- K + = K 


w. 


w 


0 


z* = 




It - 


(cos w^t - 1) 


w. 


(w n t - sin w n t) 9 
K + — - = K 


w 


0 


A direct evaluation shows that the skew-symmetric component of 
the Z matrix will remain equal to zero if the rates about all 
three axes have the same amplitude. Without the analytic solu- 
tion, the results of a simulation run or the computable solution 
with equal rate inputs might be misleading. The two's complement 
wordlength error will generally have both a symmetric and a skew- 
symmetric component. 

For a sign and magnitude number system, exact analytical 
solutions are quite difficult to obtain because of the sgn 
function. However, approximate solutions are easily obtained. 


S/M - SPSA 


Let : 


sgn (A$) = B = B . 


Then : 


• k T 
z = J B ttB 


and if Q and B commute 


Z 



dx . 
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For constant angular rates : 


Z 



S/M - SPLA 
Let : 


sgn ($) = B . 


Then : 


• k T 

z = | b t b 


and 

Z = | It 


for all angular rate inputs. 


Sign and magnitude arithmetic with a long accumulator has 
desirable error properties since the error propagation is sym- 
metric and therefore subject to compensation by an orthogonality 
correction. Short accumulator addition, on the other hand, 
results primarily in a skew-symmetric error which is undesirable. 
In both cases, there will be second-order terms not included in 
the approximate solution. 

Computable solution .- A computable solution for the direc- 
tion cosine error may be achieved by integrating the equation set 


B = fiB (169) 

E = ftE + A 

or the equation set 

B = ftB 

. T 

Z = BA . 
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The details of a computer program, which will carry out this 
integration, have been explained previously. 

Simulation results .- The validity of the statistical model 
of wordlength error has been upheld by the reports contained in 
references 6 through 10. However, it is wise to do at least a 
limited amount of simulation for a general class of applications. 

A simple validation procedure was applied to the direction cosine 
problem. A FORTRAN simulation program which numerically inte- 
grated the direction cosines using a second-order Taylor's 
algorithm was written. The simulation was done in both single 
precision and double precision. The single and double precision 
results were differenced to determine the single precision word- 
length error. The program was run on a one's complement machine 
(Univac 1108) and on a sign and magnitude machine (IBM 7094) . 

The results were than compared to the error which was predicted 
by a computable solution. Figures 21 and 22 contain plots of 
three elements of the Z matrix obtained from the simulation and 
computed solution. The sign and magnitude computer has a long 
accumulator floating point unit (although it uses a short accumu- 
lator mechanization for fixed point addition) and therefore 
exhibits predominantly symmetric error terms in the Z matrix. 

The one ' s complement machine has a predominantly skew- 
symmetric error characteristic of short accumulator addition. 

Both machines also have smaller errors of the opposite type. 

The secondary errors are less accurately estimated by the com- 
putable solution. This is because the computable solution models 
only the principal source of the wordlength error — the addition 
of the state variable increment to the state variable. If this 
error is reduced by double precision or extended precision addi- 
tion, or if a more accurate evaluation of the secondary errors 
is desired, the errors in the increment must also be included. 

This may be done by completely analogous techniques. The terminal 
values of the plots are given in Table XVI. The asterisks mark 
the errors which were predicted by the approximate analytical 
solutions . 

The simulations were also run with a continuous orthogonality 
correction. The results are given in Figure 23. They illustrate 
the suppression of the symmetric error component by the ortho- 
gonality correction. It is important to remember that the word- 
length errors are random variables. The computable solution of 
Figures 21 and 22 reflects only the expected (or mean) value of 
the error. It provides an estimate of the average error to be 
expected from a number of experiments. A simulation run is just 
one such experiment. The actual wordlength error determined by 
a simulation run will have some variation about the mean. This 
variation is measured by the covariance matrix P (t) . For many 
applications, the variations are small compared to the mean value. 
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TABLE XVI 


TERMINAL VALUES OF WORDLENGTH ERROR 


Addition 
SPSA 
Z x 
Zxx 
Zxy 

* 

Simulation 

0.605 x 10 -4 
-0.250 x 10“ 6 
0.937 x 10“ 7 

Computable Solution 

0.591 x 10 -4 
0.876 x 10~ 6 
-0.275 x 10 -6 

SPLA 




Zx 


0.181 x 10* 5 

-0.264 x 10 -6 

Zxx 

* 

0.114 x 10" 3 

0.116 x 10 -3 

Zxy 

* 

-0.177 x io~ 5 

-0.169 x 10 -5 


Under these conditions, each simulation run will closely approxi- 
mate the mean value predicted by the computable solution. This 
is the case of the simulation results of Figures 21 and 22. For 
the secondary error and when the mean value of the symmetric error 
is suppressed by an orthogonality correction, the results are 
essentially random, indicating that the variations about the mean 
are no longer negligible compared to the mean. 

The computable solution provides an estimate of the covari- 
ance matrix P (t) as well as an estimate of the mean value. For 
those applications where the standard deviation of the error is 
much smaller than the mean value, a reduction in the overall error 
can be achieved by using a symmetrical rounding procedure. The 
rounding procedure described previously required the addition of 
a small correction factor to each word before it is shortened. 
Additional machine time is required to round off numbers in this 
fashion. Reference 13 contains another rounding procedure. It 
consists of forcing the last bit of each word to be "1". This 
procedure results in a zero mean wordlength error but with a 
distribution range twice that shown in Figure 19. The simulations 
of Figures 21 and 22 did not employ a rounding operation. An 
example of symmetric rounding (applied to a different problem) is 
contained in reference 9. 




IV. SYSTEM PERFORMANCE EVALUATION 


Up to this point, the analysis of the direction cosine 
error has been concerned with the error matrices E and Z . Al- 
though the solutions provide insight into the nature of the direc 
tion cosine error and, to a limited extent, allow comparisons of 
the different methods of calculating the cosines, a more definite 
solution is desirable for evaluating system performance. 

Reference Trajectory 

A reference trajectory is specified as a time history of 
vehicle angular rates and accelerations, viz.: 

a c = vehicle acceleration in body coordinates, 

a N = vehicle acceleration in navigational coordinates, 

v = vehicle velocity in navigational coordinates, 

r = vehicle position in navigational coordinates, 

g(r) = gravitational acceleration computed from the current 
estimate of vehicle position, 

B = true direction cosine matrix, 

/\ 

B = on-board estimate of the cosine matrix. 


Except for the rather restricted case of a self -commutative 
rate matrix, a solution for the true direction cosine matrix will 
in all probability, not be known. Under these circumstances, the 
true cosine matrix B is really a fictitious quantity. However, 
it is usually possible to determine a fairly accurate matrix by 
using a sophisticated integration technique such as the fourth- 
order Runge-kutta method. This matrix will be donoted by B* . 

The nominal mission is "flown" by integrating the equations 


B* = Q* (170) 

b 

v = a N + g (r) 

r - v 
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The result is a time history of nominal vehicle position and 
velocity. If the acceleration time history is available in body 
coordinates rather than navigational coordinates, the transforma- 
tion 


- E * Ta B (171> 


may be used. The error in the cosine matrix is of no consequence 
since any actual mission will not follow the nominal trajectory 
exactly. All that is necessary is a profile reasonably close to 
the nominal to use as a reference. 

The effect of direction cosine errors can be evaluated by 
simultaneously integrating the equations 

6v = Z T a M + 6g(r) (172) 

N J 

6r = <5v 


where 6g(r) is determined from 6r and the nominal r vector. 

The Z matrix can be found (for a general angular rate input) 
(1) by finding an analytic W matrix. Then: 

R = B* T WB* (173) 


where the errors in the B* matrix generate second order errors in 
the R matrix. The Z matrix is then found by integrating 


Z = RZ + R , 


or (2) by using a computable solution for the Z matrix. 

For either method, it is necessary to make a number of runs 
for the different types of errors and algorithms. It is not 
necessary to make runs for different computational rates or word- 
lengths since these appear as parameters in the error expressions. 
The algorithm error and the mean value of the wordlength error 
may be treated as deterministic. The variation of the wordlength 
error about the mean must be given a one-, two-, or three-sigma 
value depending upon the corresponding choice for the statistical 
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instrument errors. The validity of the final results is highly 
dependent upon the choice of a realistic vehicle environment. 

For example, the rate matrix fi(t) may include sinusoidal terms 
to represent vehicle limit cycling and angular vibration. Pseudo- 
random numbers can also be generated and passed through a shaping 
filter to model random vibrations. 

The output of the error analysis program is the velocity, 
position, and attitude errors of the vehicle. The number of 
error terms can be reduced by calculating the magnitude of the 
velocity error from the vector components and using this as a 
single measure of system performance, or a figure of merit (FOM) 
may be based upon the midcourse correction velocity required of a 
spacecraft as a function of all error terms . The important point 
is that the same criteria be used for the software and hardware 
errors. The next section will outline a software design technique 
which will optimally allocate computing resources on the basis of 
the error analysis described here. 

Software Design Technique 

As the previous sections have indicated, it is not always a 
simple matter to determine the "best" direction cosine algorithm. 
The choice will be highly dependent upon the accuracy requirement, 
operating environment, and the computing resources available. 

These parameters influence the design of the entire software 
package, so the technique to be outlined may be applied to other 
programming areas as well as to the direction cosines. 

A basic requirement of the software is that it maintain an 
acceptable level of computational error in a given environment. 

The accuracy specification is often taken as a fraction of the 
hardware error for the same mission. For this reason, the soft- 
ware error should be expressed in terms of the same figure of 
merit (FOM) as the hardware error. 

The error for a software routine will have the form: 


FOM = EAL + EWL 

where 

EAL = algorithm error 
EWL = wordlength error 

FOM = figure of merit which is to be minimized. 
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The individual terms have the form: 


EWL = f*f 
1 c 

EAL = k*/f n + k*/f m 
n' c nr c 

where the k^ ' s are mission-dependent constants. For example, a 
second-order method would have one error term of the form: 

EAL = k*/f 4 . 

A C 

For a given computational frequency, the fourth-order error 
will normally be much less than the second-order error. This is 
not a meaningful comparison, however, since a second-order 
algorithm is simpler and, for an equal amount of computer time, 
it could be executed at a higher computational frequency. This 
discrepancy may be resolved by comparing the error terms on the 
basis of an equal percentage of computer time. Let ' T^ be the 
time required to make one pass through the algorithm. Then: 


is the fractional part of the computer time which must be devoted 
to the algorithm. The FOM may then be re-written as 


FOM 


k l f t 


+ k /f J 
n ' t 


k /f 
nr 


m 

t ’ 


Each algorithm will also occupy a fractional part of the 
computer memory. This may be denoted by f m . The total cost of 
implementing any algorithm is then given by the complex number: 


cost = f, + jf 
t J m 

The computational resources that have been allocated are 

F, = \ f . v and F = ^ f ^ 

t / j tK m / j mK 

K K 
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where the summations are over all the algorithms currently in- 
cluded in the software package. The resources available are 


RES = (1 - F. ) + j (1 - F*) . 

t m 


The designer may now select the algorithms which minimize the 
error at a specified cost or he may minimize the cost required to 
meet an error specification. In the latter case, a possible per- 
formance index is the complex variable 


J 



j 


1 


F 

m 


F 


m 


or the real variable 


J 




In either case, the search to minimize the performance index 
will be over a discrete set of algorithms. In programming the 
algorithms, it will often be possible to trade-off memory re- 
quirements against computational time. Thus, one basic algorithm 
may have many variations. Each variation is included as a 
separate entity. As the software package nears completion, it 
may happen that almost all the available computer memory is used 
and F m is almost unity. The performance index will then greatly 
favor routines which require little computer memory. Despite this, 
the total memory required may eventually exceed that which is 
available. In this case, another memory module must be added to 
the system. When this is done, the performance index may then 
weigh time efficiency much more heavily than memory efficiency, 
and a completely different choice of algorithms would be optimum. 
The choice of a "best" algorithm is thus complicated by the 
circumstances of its employment. In addition, the overall 
optimization will be constrained by the current extent of the 
software library. 


Electronics Research Center 

National Aeronautics and Space Administration 
Cambridge, Massachusetts, September 1968 
125-23-02-09 
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APPENDIX A 


INITIAL CONDITION 


The direction cosines propagate as 


B = fiB . 


The initial condition may always be taken as the unit matrix 
without loss of generality since, for an arbitrary initial condi- 
tion : 


B (0) = B q . 

A new matrix may be defined as 


U = BB 


T 


■ • fji 

u = bbJ = «bbJ 


U = S2U U ( 0 ) = I 


so that the solution for B and U differ only by a constant matrix 


B 


0 


APPENDIX B 


EXACT SOLUTIONS 


For a self-commutative (but otherwise arbitrary) rate 
matrix, an approximate R matrix is easily found by using the 
narrow-band approximation. For the approximation to be valid, 
it is necessary that wqIi be a small number. Since this condition 
will be satisfied for a practical application, there is not much 
incentive for finding an exact R matrix. However, this can 
always be done. From Eq. (53) of the main body of text: 


$ = i + JL. K + 

w. 


(1 - c) 


w. 


where 


s = sin WqH 
c = cos WqH. 


If a first-order Taylor series algorithm is used: 


$ = I + HK 

R = f ($ T 8 - I) 
c 

or 

R = f Fh — - (1 - c) h]k + l 1 ~ C - — 

= L W 0 J V „2 «0 

This is an exact expression for the R matrix. To find a closed- 
form solution for the Z matrix, it is necessary to integrate the 
R matrix. The integration is usually straightforward for the 
approximate R matrix determined by the narrow-band method; how- 
ever, for the exact R matrix, it may be considerably more diffi- 
cult. An exception is the constant rate case when 



H = h . 
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The resulting solution for a and B are given in the Table for 
both the first- and second-order Taylor series. The table also 
contains the narrow-band approximate solution. The exact solu- 
tion reduces to the approximate solution for small values of wgh. 
The exact solution requires the subtraction of nearly equal num- 
bers and for this reason, the approximate solutions are better 
suited for numerical computation. 

For a sinusoidal input, the integration of the exact R 
matrix becomes more difficult since it involves integration of 
terms of the form 


sin 

(m 

sin 

wt) 


cos 

(m 

sin 

wt) 

• 

The integration 

can be 

carried out by noting that 

cos 

(m 

sin 

wt) 

= J Q (m) + 2J 2 (m)cos 2wt + ... 

sin 

(m 

sin 

wt) 

= 2J-, (m) sin wt + 2J_. (m) sin 3wt 


This allows a series solution in terms of Bessel functions which 
is not restricted by the narrow-band criteria. The solution has 
little practical importance, however. 


CONSTANT RATE CASE - TAYLOR SERIES ALGORITHMS 


Order 

Exact Solution 

Narrow-Band 

Solution 

a 

3 

a 

B ! 

1 

w ( c s \ t 

2A- C s\ 

3 

" W ° t 

2 

W ° t 

M c » 0 j t 

ol 2 w n r 

Vw o 0/ 

3f 2 
c 

2f c 

2 

(_ s __ ^sV 

2/l-c s h c\ 

3 

^0 t 

4 

W ° t 

w 0^ c kw Q w 0 2 p 

°Vhw 2 "0 2 ) 

6f 2 
c 

8f 3 
c 
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APPENDIX C 


EQUIVALENT INSTRUMENT ERRORS 


Since the inertial system's accelerometers will have scale 
factor and skew (misalignment) errors and gyroscope errors may be 
expected to cause a drift of the navigational reference frame, it 
would be helpful if the computational (direction cosine) errors 
could be related to equivalent instrument errors. One reason for 
this approach is that the effects of the instrument errors on the 
position and velocity must be determined anyway, and so they may 
be presumed to be already known. A second reason is that, from a 
system standpoint, the computational error must be balanced against 
system hardware performance. Therefore, it might be convenient 
to describe the performance of the computer with the same indexes 
that are used for the inertial instruments without reference to a 
specific mission. Unfortunately, the analogy between cosine 
errors and instrument errors is somewhat tenuous. Let 


SF = the diagonal matrix of accelerometer scale factor 
errors , 

SK = the symmetric matrix of accelerometer misalignment 
angles , and 

D = the matrix representing the integral of the system 
gyroscope drift rates 


Then: 

1. Since the accelerometer misalignments result from 
construction inaccuracies , the SK matrix would be 
expected to be neither symmetrix nor skew-symmetric 
and to contain six independent error angles. How- 
ever, any matrix may be resolved into symmetric and 
skew-symmetric components. The skew-symmetric part 
of the SK matrix may be compensated by the initial 
alignment of the inertial system leaving SK symmetric 
with only three independent angles. 

2. The integral of the gyro drift rate is used for D 
since the analysis of the direction cosine matrix 
led to the definition of Z and E matrices which are 
angular error matrices. 
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The acceleration error due to the instrument errors defined 
above will be 


a^ = (I + D + SI)ag 


where 


SI = SF + SK 
<Sa^ = B T (D + S)a 3 
S'a N = B T (D + S)Ba N . 
The computation error is given by 




f 


so that an equivalence is given by 

T T 

D + S = BZ B . 


The above equation represents a congruence transformation of the 
Z matrix which results from the fact that the instrument error 
is defined in the body coordinate frame while the Z matrix 
(representing the computational error) is defined in the naviga- 
tional coordinate frame. In general, the computational and 
instrument errors are related by the time varying B matrix. An 
equivalence can now be defined by 


T T 

D = B ( -Z - Z Z)B 

ss 

T T 

S = B ( Z + Z Z)B . 
s 

When an orthogonality correction is used, the skew error will be 


97 



zero. When the angular rate matrix is self-commutative : 

T 

D = -Z - Z Z 
ss 

S = Z + Z T Z . 
s 

The D matrix can be approximated by 


which is skew symmetric. 

It is common practice to assume a constant gyro drift term 
(D will then be a linear function of time) and a constant S matrix. 
However, Z may be a more complex function of time. Therefore, any 
equivalence between instrument and computational errors can only 
be meaningful in the context of equivalent instrument coefficients 
which are allowed to be functions of time. 
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APPENDIX D 


EFFICIENT COMPUTATION 


The algorithms used in this report were written in matrix 
form which is convenient for analysis purposes. The programming 
of the equations, on the other hand, is often considerably 
simplified if it is not done in matrix form. 

With the skew-symmetric matrix of gyroscope inputs 0 [see 
Eq. (62)] a three-dimensional vector ij; may be associated where ip 
contains the three independent elements of the matrix 0. From 
Eq. (95) the exact transition matrix is seen to have the form 

$ = I + a0 + b0 2 . (D-l) 


This can be written as 


$ = I + 0 A + 0 0 0 (D-2 ) 

where 

<P A = aip (D-3) 

= bijj 


with the modification of ip applied to a vector rather than to a 
matrix . 

Approximate transition matrices have the same general form. 
From Eqs . (66) and (95), it is evident that a fourth-order 

approximation has the form of Eq. (D-2) provided that 

a = 1 - <J> 2 /6 (D-4 ) 

b = 1/2 - <p 2 /24 


where 


is defined by Eq. 


(95) . 
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The same approach can be applied to the compensated solu- 
tion of Eq. (120). Two gyroscope samples are obtained (\p lf ^ 2 ) 
for each calculation of a compensated cosine matrix. Eq. (120) 
can then be approximated by 


rp = ip 1 + ip 2 + y(ip 1 x ip 2 ) . 


(D-5 ) 


The sequence of calculations is then Eqs. (D-5), (D-4) , (D-3) , 

and (D-2) . In addition Eq. (D-2) need not be programmed in matrix 
form. The result is a fourth-order algorithm with a commutativity 
correction incorporated in Eq. (D-5) . 
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DIRECTION COSINE COMPUTATIONAL ERROR 


By John W. Jordan 
Electronics Research Center 
Cambridge, Massachusetts 


ABSTRACT 

Strapdown inertial systems replace the gimbal structure of 
more conventional inertial systems by a direction cosine matrix 
which is contained in the system computer. Since the cosine 
matrix must be updated many times a second, a detailed study of 
the computational error introduced by the on-board computer is 
necessary to determine system performance. This report is con- 
cerned with the techniques which may be employed to estimate the 
computational error and its effects on system performance. 
Although focused upon a particular problem of practical impor- 
tance, many of the techniques have a wider applicability to aero- 
space software in general. The computer-generated error is 
divided into two categories: (1) Algorithm Error which is caused 

by the approximate nature of the numerical formulas used, and 
(2) Wordlength Error which is caused by a finite computer word- 
length . 

Wordlength Error is treated from the statistical viewpoint 
of Henrici. The theory is developed for linear matrix differ- 
ential equations. The direction cosine problem is then used as 
a specific example. The results of Henrici are extended to in- 
clude different computer number systems (i.e., sign-magnitude or 
complement representations) and different organizations of the 
computer's arithmetic unit. The emphasis is upon a computable 
solution which may be applied to problems for which analytical 
solutions are either not possible or are not practical. The 
solution includes both fixed- and floating-point machines. The 
computable solution is verified by simulation. 

A differential equation is derived for the propagation of 
the Algorithm Error . Since the direction cosine matrix should be 
orthogonal, it is possible to include a correction routine in the 
airborne computer so that the cosine matrix will remain orthogonal 
despite computational error. A second differential equation is 
derived for the propagation of the Algorithm Error when an 
orthogonality correction is used. Closed-form analytical solu- 
tions are obtained for both error differential equations provided 
that the vehicle angular rate matrix is self-commutative. 

In the general case, the error is determined by a computable 
solution or by simulation. The report also considers compensation 
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techniques, the evaluation of system performance, and the "optimum" 
design of the aerospace software. All analytical solutions are 
verified by simulation. 


STAR Categories 08, 19, and 21 
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